Local Fourier Analysis of the Complex Shifted Laplacian 
preconditioner for Helmholtz problems 



Siegfried Cools, Wim Vanroose* 

Department of Mathematics and Computer Science, University of Antwerp, Middelheimlaan 1, 2020 Antwerp, Belgium 



SUMMARY 

In this paper we solve the Helmholtz equation with multigrid preconditioned Krylov subspace methods. 
The class of Shifted Laplacian preconditioners are known to significantly speed-up Krylov convergence. 
However, these preconditioners have a parameter /3 G M, a measure of the complex shift. Due to 
contradictory requirements for the multigrid and Krylov convergence, the choice of this shift parameter 
can be a bottleneck in applying the method. In this paper, we propose a wavenumber-dependent minimal 
complex shift parameter which is predicted by a rigorous /c-grid Local Fourier Analysis (LFA) of the 
multigrid scheme. We claim that, given any (regionally constant) wavenumber, this minimal complex shift 
parameter provides the reader with a parameter choice that leads to efficient Krylov convergence. Numerical 
experiments in one and two spatial dimensions validate the theoretical results. It appears that the proposed 
complex shift is both the minimal requirement for a multigrid V-cycle to converge, as well as being near- 
optimal in terms of Krylov iteration count. 

KEY WORDS: Helmholtz equation, indefinite systems, high wavenumber, Krylov method. Shifted 
Laplacian preconditioner. Local Fourier Analysis 



1. INTRODUCTION 

The propagation of waves through a material can be mathematically modelled by the Helmholtz 
equation 

on a d-dimensional subdomain ^7 c M'^. Here k{x) is the space dependent wave number that models 
change in material properties as a function of x ^Vt. For high wavenumbers k, the sparse system 
of linear equations that results from the discretization of this PDE is distinctly indefinite, causing 
most of the classic iterative methods to perform poorly. Incited by its efficiency on positive definite 
problems, the more advanced multigrid method has repeatedly received attention as possible solver 
for these systems. However, as stated in [1] and [2], a direct application of the multigrid method to 
the discretized systems will inevitably result in divergence, as both the smoother and the coarse grid 
correction tends to introduce growing components in the error 

The aim of this paper is to understand the required modifications to the multigrid method such that 
it effectively solves the indefinite linear systems resulting from the discretization of the Helmholtz 
equation. 

Over the past few years, many different methods have been proposed to solve this non-trivial 
problem, an overview of which can be found in [3]. Krylov subspace methods like GMRES 
[4] or BiCGStab [5] are known for their ability to definitely converge to the solution. However, 
Krylov methods are generally not competitive without a good preconditioner. Consequently, in 
recent literature a variety of specifically Helmholtz-tailored preconditioners have been proposed to 
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speed up Krylov convergence. In this paper we focus on the class of so-called Shifted Laplacian 
preconditioners, which were introduced in [6] and further analyzed in [7] and [8], where they 
are shown to be very efficient Krylov method preconditioners, significantly reducing the number 
of Krylov iterations on both academic and realistic problems. Furthermore, unlike the original 
Helmholtz system, the Shifted Laplacian preconditioning system can be solved effectively by a 
multigrid method. 

Similar to many other Helmholtz solvers like e.g. [9] however, the Shifted Laplacian method 
introduces a parameter which is intrinsic to the method itself, namely the magnitude of the (real 
and/or imaginary) shift. Ideally no shift would be introduced, as preconditioning the original system 
with its exact inverse causes the Krylov method to immediately converge to the solution without 
iterating. Unfortunately, the use of a (complex) shift is vital to guarantee multigrid convergence, as 
stated in [3], [6], [7] and [10]. Summarizing, one could say that the Krylov method prefers the shift 
to be as small as possible, whereas the multigrid solver only converges for a sufficiently large shift. 
This contradiction makes the choice of an optimal shift a non-trivial yet essential task in fine-tuning 
the Shifted Laplacian method. 

Note that, whereas in this paper the shift is implemented directly into the discretization matrix, 
it can also be represented by using complex valued grid distances in the discretization [11]. 
Furthermore, a complex shift can be taken into account in algebraic multilevel method, as elaborated 
in [12]. The Krylov convergence rate has been explored as a function of a complex shift in [13]. 

The current paper aims at gaining more insight in the value of the complex shift, up to current 
date, this shift is determined primeraly on a heuristical case-by-case basis, as very limited theoretical 
founding is known. Furthermore, in recent literature a discussion has risen on the question whether 
the choice of this complex shift can be independent of the wavenumber k. We will introduce the 
notion of a minimal complex shift, which we define as the smallest possible shift for multigrid to 
converge, and show that this shift must indeed depend directly on the wavenumber. Additionally, 
it will be verified that the proposed complex shift is optimal in terms of preconditioned Krylov 
iterations when exactly solving the preconditioning problem. When applying a limited amount of 
multigrid iterations to the preconditioning problem as is common practice, only approximately 
solving the preconditioning system, it will appear that our definition provides an near-optimal 
value for the complex shift parameter w.r.t. Krylov convergence. Consequently, following the given 
definition, the reader is provided a near-optimal and generally safe choice for the complex shift. 

The main tool used in this paper to analyze multigrid convergence and determine a realistic value 
for the complex shift is the Local Fourier Analysis, originally introduced by Brandt in 1977 [14]. 
Our analysis is primarily based on the LFA setting introduced in [15], [16] and [17]. Moreover, the 
analysis in this text can be seen as an extension of the spectral analysis performed in [18] and is 
even more closely related to recent work like [19] and [20], the latter combining Shifted Laplacian 
with the new technique of multigrid deflation for improved convergence. 

The results in this paper are limited to problems with homogeneous Dirichlet boundary 
conditions. Realistic Helmholtz problems often have absorbing boundary conditions implemented 
with an absorbing layer such as PML [21]. Note that these absorbing boundary conditions often 
lead to a more favorable spectrum for Krylov convergence. The analysis performed here can thus be 
considered a worst-case convergence scenario for multigrid preconditioned Krylov methods. 

This paper is organized as follows. In Section 2, we present a rigorous LFA analysis of the 
problem, ultimately allowing us to define and effectively calculate the minimal complex shift 
parameter. In Section 3, numerical results for actually solving the Helmholtz model problem with 
a variety of different wavenumbers and grid sizes are presented, which serve as a validation of the 
definition given in Section 2. These experiments verify that the theoretical complex shift is indeed 
minimal w.r.t. multigrid convergence, and we furthermore observe near-optimality w.r.t. Krylov 
convergence. Additionally, we briefly discuss the similarity between our results and the observations 
made in [22]. We conclude the paper with a short review of the conclusions in Section 4. 



Note that throughout the text, we will use the symbol l to denote the imaginary unity or to 
avoid confusion with the index designator i. 
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2. LOCAL FOURIER ANALYSIS 



2.1. Problem statement and notation 

The general aim of our research is to solve the d-dimensional indefinite Helmholtz equation on an 
open bounded domain 1] c 

— au = f onQ, (1) 

where a G Mq ^ distinctly negative constant. Observe that we denote the squared wavenumber as 
—a = k'^.^ We do not attend to boundary conditions, as Local Fourier Analysis does not take into 
account the domain boundaries, see Section 2.2. Using multigrid on a Complex Shifted Laplacian 
preconditioner, we consider the problem of iteratively solving the related (i-dimensional complex 
shifted Helmholtz equation 

— A'u -\- au = g on 1^, (2) 

where a = cr(l + I3l) G C, with complex shift parameter (3 G M+. Note that the CSL preconditioner 
was originally introduced in [6] with a somewhat more general a = a{a -\- Pl), but following the 
observations in [7] we will choose to permanently set a = 1. This is in a way a natural choice for 
a, as preserving the real shift keeps the preconditioner very close to the original problem. Equation 
(2) is typically discretized using a finite difference scheme, yielding a linear system Au = g. In 
our multigrid analysis, we assume a standard second-order finite difference discretization on an 
equidistant mesh with N gridpoints, where typically N = 2^ for some p G No and with mesh width 
h = 1/N. This renders for d = 1 a discretization matrix A^^ with stencil representation 

^i^ = ^( -1 2 + -1 ), (3) 

or, for (i = 2, a discretization matrix A'^^ given by 



12^ = ^1 -1 4 + -1 I, (4) 



which can easily be generalized for higher dimensions. Throughout the analysis, we consider the 
propagatioi 
as follows 



propagation of the initial fine grid error e^^^ through a k-grid analysis error propagation matrix 



^(m+l) ^^,^(m)^ m>0, (5) 

where designates the iteration matrix on grid level / corresponding to the multigrid cycle that 
employs k grid levels. Note that we designate the finest grid by the index / = 1, the one-level coarser 
grid by / = 2, and so on, implying a total of 2^~^+^ grid points on the l-th grid. For the well-known 
two-grid cycle, can be expressed as (see [15]) 

m2 = s'(Hh - 4a^'i!a^)s';\ (6) 

which is generalized by the following definition of the error propagation matrix Mj^ for a k-grid 
analysis (see [16], [17], [19]) 

M,' = Spill - - Mt^,)Al-_^\ll^'Ai)Sp, / = 1, . . . , - 1, 

Ml: = 0, (7) 

where Si is the smoothing operator, Ai is the discretization matrix representation and Ii is the 
identity matrix on the l-th coarsest grid Gi with mesh size hi = 2^~^h {I = 1, . . . ,k). For the one- 
dimensional problem, the l-th coarsest grid Gi is given by 

Gi = {xen\x = Xj= jhu j G Z}, (8) 



'''To avoid unnecessary terminological complications, in this text we loosely refer to a as 'the wavenumber'. However, 
the reader should keep in mind that we hereby intrinsically designate the negatively signed squared wavenumber — /c^. 
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whereas in 2D it is defined as 



= {x G I x = {xj^.Xj^) = {jihij2hi), jij2 e (9) 

Furthermore, : Gi Gz+i and : Gi are the restriction and prolongation 

operators, respectively. In this paper, we will use full weighting restriction and linear interpolation 
as the standard intergrid operators. 

2.2. Basic principles of ID Local Fourier Analysis 

We briefly sketch the key ideas behind Local Fourier Analysis. LFA is based on the assumption that 
both relaxation and two-grid correction are local processes, in which each unknown is updated using 
only the information in neighbouring points. Furthermore, boundary conditions are neglected by 
extending all multigrid components to an infinite grid. It is presumed that the error ej^^ on the l-th 
coarsest grid can be written as a formal linear combination of the Fourier modes (pi{0, x) = e'^^^l^^ 
with X ^ Gi and Fourier frequencies 6> G M. These frequencies may be restricted to the interval 
6 = (— 7r,7r]cMasa consequence of the fact that for x e Gi 

^i{0^27r,x)=^i{e,x). (10) 

The set of l-th grid Fourier modes is typically denoted 

£i = span{ipi{0,x) = e'^^^^^\xe Gi, e 9}. (11) 

The Fourier modes are known [23] to be formal eigenfunctions of the operator Ai. More precisely, 
the general relation Ai(pi{0,x) = Ai{0)(pi{0,x) holds, where the formal eigenvalue Ai{0) is called 
the Fourier symbol of the operator Ai. Given a so-called low Fourier frequency 6^ G (— 7r/2, 7r/2], 
its complementary frequency 0^ is defined as 

=0^ - sign(i9^)7r (12) 

Using this notation, the following important property can easily be derived (for 1 < I < k) 

^i{e\x) = ^i+r{2e\x) = ^i{0\x), X e Gi^i. (13) 

The Fourier modes (pi (6>^, x) and (pi {0^ , x) are called (/ + l)-th level harmonic modes. These Fourier 
modes coincide on the (/ + l)-th coarsest grid, where they are represented by a single coarse 
grid mode with double frequency (pi^i{20^ ,x). In this way, each low-frequency mode (pi{0^,x) 
is naturally associated with a high-frequency mode (pi{0^,x) on the l-ih grid. It is convenient to 
denote the two-dimensional subspace of Si spanned by these (/ + l)-th level harmonics as 

Sf = spm{^i{0', •), •)}, l<l<k, (14) 

with 0^ as defined by (12). Note that every coarse-level subspace Si^-^ can be decomposed into 
spaces spanned by finer- level harmonics, as (13) implies 

^f+i = ^f^^ U 1 < / < A: - 1. (15) 

The significance of these spaces f f is that they are invariant under both smoothing operators and 
correction schemes under general assumptions. Due to the previous observations, throughout the 
Fourier symbol calculation one can assume without loss of generality that each l-th grid error ej^^ 
can be decomposed into components e[^^ = e[^^ {xj ) that consist of a single Fourier mode (pi (6>, Xj ) 
and thus can be represented as 

= A(^) ipi{0, Xj) = A(^) e^^'^ 6> G e, j G Z, m > 0, (16) 
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where the amphtude A^^^ changes as a function of the iteration m. Moreover, using this expression 
together with the error relation in (5) and the stencil representation of the operators forming the error 
propagation matrix , it can be derived that the error amplitudes in two subsequent iterations are 
related by an amplification factor function Qi{0^ a, (3) 

^(m+i) ^ gi(^o,cT,p)A^'^\ m > 0, (17) 

which describes the evolution of the amplitude A*^^^ through consecutive iterations. It is our aim to 
elaborate an analytical expression for the amplification factor on the finest grid Gi{0, a, (3), which 
can be considered a continuation of the spectrum of (see [15], [10]). For notational convenience, 
we tend to denote Gi{0, a, (3) as G{0, a, (3), dropping the fine grid index. This amplification factor 
can be computed by calculating the Fourier symbols of each component of the Gi-grid error 
propagation matrix . After the symbols are structured in so-called eigenmatrices for each 
component, combining harmonic frequencies, and the appropriate products of these matrices are 
taken, a 2^~^ x 2^~^ eigenmatrix representation of the ID error propagation matrix is 
obtained. The amplification factor function Gi{0, a, (3) can then readily be computed as the spectral 
radius of (see [16], [23]) 

G{0,a,f3) = p{M^) = max|A(M^)| (18) 

A well-known condition for a general iterative methode with given iteration matrix, say M, to have 
convergence, is that p{M) < 1. This convergence condition can be generalised to the following 
demand on the Gi-grid amplification factor 

max G{0,(j,p) < 1, (19) 

clearly representing the analogue condition for a multigrid cycle. We denote the amplication factor 
as a function of the Fourier frequency 6, the wavenumber cr and the complex shift parameter /3. 

Definition of the minimal shift. With the help of the amplification factor we define the minimal 
complex shift parameter /3min(cr) as a function of the wavenumber a as 

/^min := argmin < max GiO, a, /3) < 1 > . (20) 

/3>0 L J 

This definition can be interpreted through (19) as the smallest possible complex shift required 
for the multigrid method to converge, i.e. it is the smallest value of (3 for which every single 
eigenmode of the error is reduced through consecutive multigrid iterations. Additionally, the 
numerical experiments presented in Section 3 will show that the complex shift parameter /3min as 
defined here is near- optimal for any multigrid-preconditioned Krylov method in terms of iteration 
count. This means that when the preconditioner is inverted exactly (up to discretization error) using 
a sufficiently large amount of multigrid steps, a minimal number of Krylov steps is required when 
choosing the value of the complex shift equal to /^min- 

23. The ID Fourier symbols 

In this section we will effectively calculate the Fourier symbols of the different component 
operators of a ID /c-grid scheme. 

Discretization operator. The evolution of the error under the discretization operator Ai can be 
calculated from its stencil representation (3), yielding 

e(?+^) = i f-elTL + (2 + dhf) el?) - e\^L V .76 Z. 



hi . 

Using expression (16), the amplitude relation is found to be 
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and hence the discretization operator Fourier symbol Ai (0) is 



M0)= (-;^cos^+^+a). (21) 



Restriction operator. Using an analogous stencil argument, the error propagation under 
application of the full weighting restriction operator ij^-^ can be derived to be 

(m+l) _ 1 (m) 1 (m) 1 (m) . ^ 

Substituting the error components by (16), one obtains the Fourier symbol ij^^iO) 

il^\0) = ^{cos0^1). (22) 
Interpolation operator. The linear interpolation operator propagates the coarse grid error 

as 

(m+l) _ (m) (m+l) _ 1 (m) 1 (m) • ^ 

where the first equation holds if x e Gi fl G^+i and the second if x G G/\G^+i. Once again 
substituting the error components by (16) and combining the above cases, the following expression 
for the Fourier symbol i\^i{0) is obtained 

iU{e) = \{^ose + i). (23) 

Note that full-weighting restriction and linear interpolation are dual operators, yielding the exact 
same Fourier symbol i\^^{6) = i\j^-^{6). 

Smoothing operator. As a smoother we use standard weighted Jacobi relaxation. It is easily 
derived (see e.g. [15]) that the cj-Jacobi relaxation matrix Si = Ii — uoD^^Ai relates the error in a 
gridpoint Xj G Gi through subsequent iterations by 

Presuming the error is of the form (16), one obtains the following amplitude relation and cj- Jacobi 
smoother Fourier symbol Si{0) 



S,iO)=[l-.+ ^^cos0). (24) 

Once calculated, the ID Fourier symbols of each component of can now be structured in 
eigenmatrices, representing the action of a l-ih grid component on the subspace of (/ + l)-grid 
harmonics, where / = l,...,/c — 1. This matrix-building will be performed explicitely in the next 
sections for specific values of k to compute the Gi-grid amplification function G{0, a, (3) and, from 
this function, the minimal complex shift parameter Pmin- 

2.4. A ID Local Fourier Analysis of the weighted Jacobi smoother 

To get some insight in the problem, we will initially perform a very basic Local Fourier Analysis of 
the cj-Jacobi smoother. It is well-known that the smoother, forming an essential part of any multigrid 
scheme, is unstable for the indefinite Helmholtz problem. This instability is caused by the smoothest 
eigenmodes which have negative eigenvalues and consequently diverge under the action of the 
smoother, as described to some extend in [1] and [2]. To make up for smoother instability, one could 
determine a complex shift based purely on the divergence of the eigenmodes under application of 
the fine-grid smoother operator Si . Hence, recalling that every two-dimensional subspace of coarse 
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Figure 1 . The minimal complex shift parameter /3min from the ID smoother analysis with N = 64 and 
(jj = 2/3 as a function of the wavenumbera. Solid black line: analytical lower limit (26) corresponding to 
Fourier frequency = 0. Solid grey line: analytical lower limit (27) corresponding to the frequency = n. 



grid harmonics C S with 0^ e (— 7r/2,7r/2] is left invariant under the actions of the smoother 
operator Si , the fine grid error propagation eigenmatrix can be written as 











(25) 



with Si{0) being the fine grid smoother symbol derived in (24). Since the matrix is diagonal, the 
calculation of the spectral radius in (18) is easy. This can then be substituted in definition (20) to 
obtain the minimal complex shift Pmin- Choosing this value for /3, one ensures that all eigenmodes 
converge under the action of . The result is shown in Figure 1 , where the minimal complex shift 
parameter (3min is plotted as a function of the wave number a for N = 64 gridpoints. 

One observes from Figure 1 that in order to compensate for smoother instability, values of /3min 
are generally very large with limcr^o Pmin = +oo, suggesting exteremely large complex shifts are 
needed to eliminate smoother instability. Due to the simple form of the eigenmatrix (25), the lower 
limit Prnin cau easily be calculated analytically. The frequency maximizing \Si{0)\ over (— 7r,7r] 
can be found by setting the first derivative of \Si (0) \ equal to zero, rendering 



-8cj^ cos 



(9 sin (9 - 4a;(l - cj)(2 + ah'^) sin = 0, 



which reveals = 0, = tt and = ± arccos{— (1 — uj) {2 -\- a h'^) /2uj} SiS local extrema. A second 
derivative check confirms ^ = and = ir to maximize the smoother symbol \Si{0)\. The value of 
Pmin can be derived by substituting both maxima in inequality (19). Presuming 6> = the expression 
for Q reduces to 

2u 

1 — UJ ' 



2 + dh? 



< 1, 



which can be elaborated using the definition a = a{l + (3l) to yield a lower limit on the complex 
shift 



l</3. 



(w - 2)a/i2 

A similar lower limit can be derived for the maximum in 6* = tt, which eventually comes down to 



(26) 



w(4 + (jh?f - 2(2 + cr/i2)(4 + ^/j2) 



(w-2)(cr/i2)2 



</5. 



(27) 



The minimal complex shift fiaiia is now defined as the minimal value of /? satisfying both 
inequalities. Note that from the left-hand of inequality (26), it readily follows that Xim^^Q = 
+00, as we observed from Figure 1. 
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Figure 2. Amplification factor function G{0,a, (3) for the ID two-grid analysis with N — as a function 
of 6 G (— 7r/2,7r/2] forhxed wavenumber and complex shift. Two leftmost hgures: wavenumbera = —500 
and complex shift parameter /3 = 0.02 (a) and /3 = 0.04 (b). Two rightmost figures: a — —5000 and 

/3 = 0.2 (c)and/3 = 0.8 (d). 



Although this short analysis of the smoother is useful to obtain some intuition, the value of 
the minimal complex shift /3min is clearly overestimated since we only consider the smoother 
operator, while in a multigrid setting the smoothest error components are removed by the coarse 
grid correction. Indeed, to guarantee multigrid convergence a much smaller value of (3 may be 
chosen compared to the shift suggested by Figure 1. In the next sections the second component of 
the multigrid method, the correction scheme, will be taken into account to attain a realistic curve for 
Pmin as a function of the wave number. 



2.5. A ID two- grid Local Fourier Analysis 

A more realistic curve for the minimal complex shift of a multigrid cycle is obtained by considering 
a basic 2-grid analysis with only one presmoothing step, for which the fine grid error propagation 
matrix is given by 

Mi = {h-liA-'lfA,)Si (28) 



As suggested earlier, every two-dimensional subspace of coarse grid harmonics Ef C E with 
0^ G (— 7r/2, 7r/2] is left invariant under the actions of both the smoother operator Si and the pure 
2-grid correction operator I2 
Ml on El is given by its 2 x 2 eigenmatrix 



/2^2 ^^1^1- The action of the total 2-grid error propagation matrix 



Ml 



h 



iim 



D' 



D 



(29) 



We use the superscript-^ notation to designate the operation that transforms a vector into a diagonal 
matrix by placing the entries of the superscribed vector along the main diagonal. The spectral 
radius of this expression, and thus the function ^(6>, a, can easily be calculated analytically or 
numerically. 

In Figure 2 we show the amplification factor for two choices of a. In panel (a) and (b) of the 
figure, where a = -500, we notice the appearance of a resonance that is caused by the coarse grid 
correction. Only this one resonant mode tends to diverge, whereas for the majority of the frequencies 
the amplification factor is significantly smaller than 1 . The appearance of such a resonance was 
already discussed in [3] and originates from the inversion of the coarse grid discretization symbol 

(2^^) in (29). Using expression (21), the frequency corresponding to the semi-asymptote can 
be approximated by 



± arcsin \ —ah^/^ = ± arcsin ^ —ah^ ^ 



(30) 

forcing the real part of the denominator of Q, i.e. the real part of the symbol A2 (26>^), to zero and thus 
maximizing the function Q{Q, a, /3) for fixed values of a and /3. Note that approximation (30) only 
makes sense when \a\ < This discussion suggests that for small values of |cr| the divergence 



8 



1 



0.9 



I 0.8 




-10000 -8000 -6000 -4000 -2000 
Wave number a 



Figure 3. The minimal complex shift parameter /3min from the ID two-grid analysis with N = 64 and 
(jj = 2/3 as a function of the wavenumber a. Black and grey points represent experimentaly measured 

values of /3min (see text). 



of the coarse grid correction scheme is to be the determining factor in the choice of the complex 
shift /3min. 

Alternatively, for high values of |cr| the maximum of G{0,a,(3) appears to be spread broadly 
around ^ = 0, as one perceives from Figure 2 (c) and (d), which implies a large range of smooth 
eigenmodes tends to diverge. The maximum here is not due to the coarse grid correction, but 
originates primarily from the divergence of the smoothing operator. To substantiate this statement, 
we consider the smoothing operator Fourier symbol Si{0) calculated in (24). Presuming —A/h? < 
a < —l/h^, approximation (30) no longer holds (no semi-asymptotes are generated), and it can 
easily be established using derivative arguments that without a complex shift (i.e. with a = a), 
expression (24) reaches a maximum in either = (if a > — 2 //i^) or its complementary frequency 
= TT (if a < —2/h'^). Note that the maximization of Si{0) was already discussed in the previous 
section. Additionally, one can verify that the absolute value of this maximum is always larger than 
1, implying divergence of the smoothest mode when applying the smoothing operator Si. It is clear 
that for large values of | a |, the main factor determining the value of /3min is the divergence caused 
by the smoother operator, rather than the coarse grid correction. 

From analytical expression (29) for the amplification factor ^(6>, a, one can now numerically 
calculate the value of /3min given a fixed wavenumber a. The maximalisation over e & is done 
using an equidistant discretization of the frequency domain with sufficiently small frequency- 
step. Subsequently, the minimal complex shift argument (3min is computed using 10 steps of the 
elementary bisection method, providing an accuracy of at least 3 decimals to /^min- Figure 3 
represents the minimal complex shift parameter Pmin as a function of the wavenumber aforN = 64 
gridpoints. Notice the significant increase of the curve around a = — corresponding to the 
difference in the underlying component generating /3min as discussed above: for wavenumbers 
cr < the divergence of the smoothing operator is the determining factor in selecting /3min, 

whereas for a > — the application of the coarse grid correction operator is decisive. 

Furthermore, we added to the figure some results for Pmin based on the experimentally measured 
convergence rate. The gray points in overplot represent these experimentally measured values of 
Pmin, determined by numerical calculation of the asymptotic convergence factors for the two-grid 
scheme on a random initial fine-grid error e^^\ whilst subjecting those factors to criterion (19). The 
experimental values are perfectly matched by the theoretical curve. The black points represent some 
similarly computed experimental values for a full V-cycle. 

It is clear from Figure 3 that the aforementioned analysis yields accurate results for the two-grid 
scheme. However, it is also clear that the addition of multiple coarser grids completely alters the 
choice of /3min for small values of |cr|. To obtain a sufficiently accurate simulation of the full V- 
cycle, a higher order k-grid scheme should be applied. In particular, we would like to guarantee 
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a near-to-exact theoretical prediction of /3min around a = — (0.625//i)^ for the V-cycle, respecting 
the kh < 0.625 criterion from [24] for a minimum of 10 gridpoints per wavelength on the standard 
[0, 1] -domain. From this perspective, a 4-grid LFA analysis will appear to be satisfactory (see next 
section). We emphasize however that this paper primarily focusses on the iterative solution, rather 
than the accuracy of the discretization. 



2.6. A ID 3 -grid and 4-grid Local Fourier Analysis 

In this section we extend the two-grid LFA analysis to a more general /c-grid analysis, motivated by 
the inaccurate results of a two-grid analysis for low values of |cr| . For notational purposes, we restrict 
ourselves to a rigourous explanation of the /c = 3 case, however the 4-grid analysis is completely 
analogous. Results are shown in Figure 4 for both /c = 3 and k = A. 
The 3 -grid fine grid error propagation matrix is given by (7) to be 



Ml = (/i - Il{h - {h - llA-^llA2)S2)A-^llA^)S^. 



(31) 



For any 0^ e (— 7r/2, 7r/2], the 3-grid operator leaves the 4-dimensional subspace C f of Gs 
harmonics invariant, implying its eigenmatrix Mf is a 4 x 4 square matrix 
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(32) 



where we introduce as a short notation 0^'^ = 0^/2 and 6>^'^ = 0^/2, and additionally 6>^'^ = 
0^/2 -\-7T and 0^^^ = 0^/2-\-7r. The amplification factor function Q{6^cf^P) is per definition the 
spectral radius of the eigenmatrix Mf . 

We remark that the resonances in the function Q discussed in the previous section indeed reappear 
in the three-level approximation. However, approximating the resonant frequency analogue to (30) 
becomes a significantly non-trivial matter fork > 2 due to the inversion of multiple Fourier symbols. 

Given a wavenumber a, the value of the corresponding minimal complex shift /3min can then 
be calculated numerically as described in the previous section. The theoretical results from the 
3-grid and 4-grid analysis are shown on Figure 4 for A/' = 64. The black dots again represent 
experimentally measured values of /3min for a full V-cycle. Note that the theoretical curve from 
the 4-grid analysis precisely matches the experimental values around a = — (0.625//i)^ = —1500. 
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2. 7. 5a^/c principles of 2D Local Fourier Analysis 

The Local Fourier Analysis performed in the previous sections can easily be extended to 2D 
problems. In this case the error ej^^ can be written as a formal linear combination of the 2D Fourier 
modes (pi{6^x.) = e'^^^i^^i+^^^cjs)/^' with x = (xj^^Xj^) G G/, which is now defined by (9), and a 
couple of frequencies = (^i, 6^2) G 9^ = (— tt, tt]^ C M^. In analogy to (1 1), the subspace of l-th 
grid Fourier modes is denoted 

£1 = span{(^,(x,^) = e^(^i-.i+^2a.,,)/^. 1^ ^ e = (^1,^2) G 6^}. (33) 

Considering a low frequency = (6>J^,^2^) G (-7r/2, 7r/2]^ and defining its 2D complementary 
frequencies as 



(34) 



one can derive the harmonics property for 1 < I < k 

(^K^'',x) = ^i{0'',^) = (^z+l(2^o^x) = ^i{0'\^) = ^i{0'\^)^ X e Gi^i. (35) 

The low frequency Fourier mode (pi{6^^^ x) and three high frequency modes (pi{6^^^ x), (pi{6^^^x.) 
and (pi{6^^^x.) coincide on the (/ + l)-th coarsest grid and thus are called (/ + l)-th level harmonics. 
Analogous to (14), we denote the four-dimensional subspace of £1 spanned by these (/ + l)-th level 
harmonics as 

Sr =mnWi{O'',-),M0'',-),y^i{O'\-),MO'\-)}, l<l<k. (36) 

Again, these spaces are invariant under general smoothing operators and correction schemes. Given 
these notations, the 2D Local Fourier Analysis itself is completely similar to the ID case, with every 
l-th grid error component e^^j^ ^^i^g represented as a single Fourier mode 

^S,.2)=^^"^^^(^'(^^-^^-^))' juheZ, m>0, (37) 

from which a relation between the amplitudes in subsequent iterations can be derived 

^(m+i) ^ g^(6/,a,/3)A(^\ m > 0. (38) 

The finest grid amplification factor Gi{0, a, (3) of the 2D operator can be calculated according 
to (18), where now represents the 4^~^ x 4^~^ eigenmatrix representation of . The Fourier 
symbols defining this eigenmatrix are derived in the next section. 

Definition of the minimal shift. Once calculated, the amplification factor function allows us to 
compute the minimal complex shift parameter, defined in 2D as 

Pmin '= argmin \ max g{0, a, ^) < 1 > , (39) 
/3>o {oee'^ ) 

which is the analogue of definition (20) for the two-dimensional case. 
2.8. The 2D Fourier symbols 

Below we briefly sum the Fourier symbols of the different component operators of the 2D /c-grid 
operator . The elaboration of these symbols is completely similar to the calculations in Section 
2.3, and is omitted here in favour of readability. 
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Figure 5. The minimal complex shift parameter /3min from the 2D 3-grid (left) and 4-grid (right) analysis 
with iV = 64 and uo — 2/?> as a function of the wavenumber a. 



Discretization operator. Using stencil representation (4) and expression (37), the discretization 
operator Fourier symbol Ai{0) can readily be found to be 



MO) 



COS 6'i — -TT cos O2 
hf 



(40) 



Restriction operator. Using an analogous stencil argument, the Fourier symbol of the 2D full 
weighting restriction operator ij^^ can be derived as 



i\'^^{9) = -(cos6>i cos6>2 +cos6>i +cosi92 + 1). 



(41) 



Interpolation operator. As is the case in the ID analysis, the Fourier symbol of the 2D linear 
interpolation operator can be derived to yield exactly the same expression as the full weighting 
restriction, i.e. 



ij^^ (^) = i (cos Oi cos O2 + cos 6>i + cos 6>2 + 1) . 



(42) 



Smoothing operator. Presuming the error is of the form (37), one obtains the following cj-Jacobi 
smoother Fourier symbol for the 2D case 

2uj 



Si{e) = 1 



{cosOi + cos 6^2) • 



(43) 



2.9. A 2D 3-grid and 4-grid Local Fourier Analysis 

Once calculated, the 2D Fourier symbols of each component of can again be structured in 
eigenmatrices, representing the action of the l-th grid components on the subspace of (/ + l)-th 
level harmonics, where / = l,...,/c — 1. The product of these matrices yields the eigenmatrix , 
from which the amplification function G{0,a,l3) can be calculated. For any given a, the value of the 
minimal complex shift parameter (3min can then be computed numerically as described in Section 
2.5. We have effectively performed the 2D computation for = 2, 3 and 4. The results from the 
3-grid and 4-grid analysis for N = 64 are shown by Figure 5. Clearly, the observations made in the 
ID case are again visible here. 

We notice a significant increase of the (3min value for wavenumbers a = —2/h'^ and smaller. 
For these wavenumbers, the maximum of the 2D amplification factor is determined solely by 
the divergence of the smoother operator. For values of a > — 2//i^, the resonances caused by the 
correction scheme that appear in the ID problem near a single frequency given by equation (30), 
now appears for a range of couples {9 1^62) for which the real part of the (combination of) coarse 
grid symbol(s) is approximately zero. As in the ID case, a 4-grid analysis accurately simulates /3min 
values extracted from V-cycle experiments up to wavenumbers as small as a = —1500. 
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Figure 6. The minimal complex shift parameter /3min from the 2D 4-grid analysis with N = 64 and 
cj = 2/3 as a function of the wavenumber a, where v — I (a), v — 2 (b), v — ?, (c) and v — A (d). 




Figure 7. The minimal complex shift parameter /3min from the 2D 4-grid analysis with N = 64 and v — 1 
as a function of the wavenumber a, where uo — 0.2 (a), uo — 0.4 (b) and uj — 0.6 (c). 



2.10. Extensions and general remarks 

Note that all previous results were constructed under the assumption that only one pre- or 
postsmoothing step is applied, i.e. v = ^ ^2 = However, often multiple smoothing steps are 
used to obtain a more accurate or faster converging iterative solution to the given problem. The 
minimal complex shift obviously depends on the number of smoothing steps cf. (7). This 
observation is depicted in Figure 6, which shows the /3min-curve for = 1, . . . , 4. The instability 
of the cj-Jacobi smoother operator, caused by divergence of the smoothest modes as described in 
[1] and [2], requires the complex shift to rise significantly for some wavenumbers cr when applying 
multiple smoothing steps. Altering the number of smoothing steps clearly has a significant impact 
on the /3min-curve. As a general tendency, one could state that the value of /3min rises (at most) 
linearly as a function of the number of smoothing steps applied. 

Another parameter of the analysis is the weight G [0, 1] of the Jacobi smoother. When applying 
the cj-Jacobi smoother, it is convenient to use the standard = 2/3 weight, which is known to be 
optimal for a ID Poisson problem (see [15]). However, for the more general Helmholtz equation 
with cr 7^ 0, it can be shown (see [1]) that the optimal Jacobi weight for the ID problem is given by 



2 



(44) 



implying the smoother weight value should be smaller when considering larger values of |cr|. As 
shown by Figure 7, altering the smoother weight does not have more than a marginal effect on the 
/^min -curve, causing it to rise only slightly as uj increases. -l- 

Furthermore, a comment should be made on the discretization-dependency of the theoretical 
curves. Note that the value of /3min is indeed dependend on the finest-mesh stepsize h = hi. 



*We remark that function values /3min for wavenumbers in the region \(t\ < 1000 should not be taken into account, as 
the 4-grid scheme is not guaranteed to correctly predict these values - see preceding sections. 
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However, it is clear from the Fourier symbol calculations that this value never appears separately 
from the wavenumber a, i.e. /3min intrinsically depends on the product gK^. Hence, the value of /3min 
remains unchanged as long as crh'^ is constant. In practice, this implies that the /3min -curve for N 
evaluated in a has the exact same value as the curve for 2N evaluated in 4cr. Doubling the number 
of finest-grid points therefore implies a stretching of the /3min -curve over the wavenumber domain 
by a factor of four. 

To conclude this section, we present a short discussion on the choice of /3min that is currently 
being used in the literature. A widely accepted choice for the complex shift parameter is /3 = 0.5, 
which was first introduced in [7] and has been used ever since by many researchers in the field. In 
this paper by Erlangga, Oosterlee and Vuik, one reads that 

"The preconditioner of choice in this paper is based on the parameters {Pi^h) = (1,0.5). 
(... ) For values (32 < O.b it is very difficult to define a satisfactory converging multigrid 
F(l,l)-cycle with the components at hand. They are therefore not considered. (...) From 
the results in Table 6 we conclude that the preferred methods among the choices are the 
preconditioners with Pi = 1. (. . . ) Fastest [multigrid preconditioned Krylov] convergence is 
obtained for (/3i,/32) = (1,0.5)." 

where Pi and P2 refer to the real and complex shift parameters, designated in this text by a and 
P respectively. These conclusions where drawn from a variety of numerical experiments with a 
fixed wavenumber a and mesh width satisfying kh = 0.625. The results can be compared to the 
theoretical value for /3min suggested by Figure 6(b) at a = —1600. Note that we apply a V-cycle 
and use the smoother weight = 2/3, contrary to the F-cycle and uj = 0.5 used in [7], however 
these differences only minorly affect the Prnin curve. It is clear that for all wavenumbers a > —1600 
(indicating the use of at least 10 gridpoints per wavelength), the corresponding minimal complex 
shift is indeed smaller than 0.5. In the case where a = —1600, Pmin equals 0.2, suggesting an 
even smaller complex shift may be used. Hence, respecting the kh < 0.625 criterion, the choice 
for (3 = 0.5 always guarantees V(l,l)-cycle convergence. The supposed near-optimality of /3min 
with respect to the number of Krylov iterations is discussed in the next section. 
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3. NUMERICAL RESULTS 

In this section, we present some experiments that will be used to confirm and validate the theoretical 
results obtained in Section 2, as well as provide us with some valuable insight on the definition 
choice of the minimal complex shift as stated above. 

3.1. Minimality of the complex shift w.r.t. multigrid convergence 

To validate the results from the previous section, we use a multigrid algorithm on a Complex Shifted 
Laplacian preconditioner for some well-known Krylov methods (see [7]), which we apply to the 
discretized model problem (1) with a standard all-one right-hand side / = 1. This leads to an outer 
Krylov iteration, where in each iteration a number of ji two-grid or V-cycle inner steps are applied 
to solve the preconditioning system. Two standard Krylov methods for solving this model problem 
considered here are GMRES [4] and BiCGStab [5]. The number of multigrid steps per outer iteration 
is denoted /i, where e.g. /i = 5 indicates five inner multigrid iterations are used within every Krylov 
step. As initial guess for the solution, we use a standard all-zero vector vq. 

A first result is displayed in Figure 8, where the number of Krylov iterations is plotted for different 
choices of //, the number of two-grid iterations per Krylov step. A tolerance of 10~^ is used and 
the maximum number of Krylov iterations is capped at 64, being the number of unknowns in the 
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Figure 10. V(l,0)-cycle-preconditioned GMRES iteration count with fi — 30 (colour) as a function of the 
wavenumber a and the complex shift /S for the 2D model problem with N = 32 (b), in comparison to the 
theoretical 4-grid LFA complex shift parameter /3min (a). 
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Figure 11. Two-grid-preconditioned Krylov method iteration count for a — —1000 (top) and a = —4000 
(bottom) in function of the complex shift /3 for the 2D model problem with N = 32. The applied method is 
GMRES with /i = 3 (a), /i = 5 (b) and /i = 10 (c). Minima can be found at (top / bottom) /3 = 0.36/0.52 

(a), /3 = 0.40/0.70 (b) and /3 = 0.42/0.72 (c). 



system. For most values of f3 < /3min the number of Krylov iterations required to reach the solution 
is excessively large. Note that for some (3 < Prnin the destructive effect of the divergent modes on the 
global convergence appears only when sufficiently many two-grid iterations are applied. This is due 
to the fact that divergence is a limit concept, implying the error may not increase during the initial 
mo multigrid iterations. However, for all shifts (3 < (3min the number of Krylov iterations increases 
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Figure 12. Two-grid-preconditioned Krylov method iteration-minimum- 13 as a function of the wavenumber 
a for the 2D model problem with N = 32. The applied method is GMRES with /i = 3, /i = 5 and /x = 10. 
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Figure 13. V(l,0)-preconditioned Krylov method iteration count for a = —64000 as a function of the 
complex shift ^ for the 2D model problem with N = 256. The applied method is GMRES with ji — 3 
(a), fi = 5(b) and /i = 10 (c). Minima can be found at /3 ^ 0.26 (a), /3 = 0.34 (b) and /3 = 0.34 (c). 



dramatically as /i grows larger, in which case we approach the theoretical curve. These experimental 
results confirm that choosing the minimal complex shift at least as large as /3min always ensures a 
safe choice for /3, independently of the number of multigrid cycles performed. The experiment can 
easily be extended to a full V(l,0)-cycle in ID and 2D as displayed by Figure 9, supporting the 
theoretical results shown on Figure 4 and 5 for the higher-order k-grid schemes. 

We can state that /3min is indeed minimal with respect to multigrid convergence, meaning it can 
be considered a lower limit to guarantee convergence. Consequently, for a large (infinite) number of 
multigrid applications, /3min is the smallest possible complex shift for any multigrid preconditioned 
Krylov method to converge. 



3.2. Near-optimality of the complex shift w.r.t. Krylov convergence 

Another important observation can be made, relating the minimal complex shift parameter /3min 
to Krylov convergence behaviour. For practical purposes, we consider the slightly smaller N = 32 
model problem, with theoretical 4-grid /3min -curve as shown by Figure 10(a). 

Temporarily focussing on a specific wavenumber a. Figure 11 shows the number of Krylov 
iterations as a function of the complex shift (3 for three different numbers of multigrid applications. 
Note that for {3 <C /3min the number of Krylov iterations is typically large as discussed in the 
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(a) 



Figure 14. Two sketches of the unshifted 2D discretization matrix spectrum for wavenumbers around a = 

(a) and a = -4//i? (b). 

previous section. Additionally however, one clearly observes the existence of a minimum number 
of Krylov iterations corresponding to a certain complex shift. It can be derived from the figure 
that the minimum is reached at /3 = 0.36 (a), /3 = 0.40 (b) and /3 = 0.42 (c) for a = -1000 and 
/3 = 0.52 (a), /3 = 0.70 (b) and /3 = 0.72 (c) for a = -4000. These shifts are optimal for the Krylov 
convergence, in the sence that they reduce the global number of iterations to a minimum. Note that 
the theoretical minimal complex shift for a = —1000 equals /3min ^ 0.42, whereas for a = —4000 it 
is /^min ~ 0.72. It appears that for a large number of multigrid iterations, the ' iteration-minimum- /3' 
approximates the theoretical Pmin- 

A similar tendency applies to other values of the wavenumber, as shown by Figure 12, where 
the complex shift corresponding to the minimum amount of Krylov iterations is plotted as a 
function of the wavenumber a. Comparing this 'iteration-minimum-/^' to the theoretical (3min- 
curve, one observes that the latter is being approximated by the iteration-minimum curves for 
increasing numbers of multigrid applications. Consequently, we call /3min near-optimal with respect 
to the Krylov convergence, implying that for sufficiently large numbers of multigrid iterations, 
the global number of Krylov iterations corresponding to P^in will be minimal. However, when 
approximately solving the preconditioning problem using only a small number of multigrid cycles, 
as is common practice, the iteration-minimum- /3 may be slightly smaller than /3min, hence the term 
'near-optimality'. Nonetheless, choosing P = /3min as the complex shift provides a 100% multigrid 
stable and low-Krylov iteration solution method. 

As an extra validation of these results. Figure 13 shows the number of Krylov iterations 
as a function of (3 for a more realistic 256 x 256-grid 2D model problem preconditioned by 
different amounts of V(l,0)-cycles. Again, the near-optimality of (3min is clearly visible, as the 
shifts minimizing the number of Krylov iterations /3 = 0.26 (a), /3 = 0.34 (b) and /3 = 0.34 (c) 
approximate the theoretical minimal complex shift parameter /3 ^ 0.34. 

3.3. General remark on the number of Krylov iterations 

An additional remark can be made here regarding the number of Krylov iterations. As demonstrated 
by experimental Figures 8, 9, 10(b) and 11, the general number of Krylov iterations required 
to accurately solve the 2D model problem with a fixed complex shift P highly depends on the 
wavenumber a. This can clearly be observed by looking at the ^-axis scale of Figure 11. The large 
differences in Krylov iterations are due to the indefinite nature of the problem, causing the number of 
iterations to drop suddenly for values of a around — 8//i^, i.e. where the 2D problem turns negative 
definite, as Figure 10(b) illustrates. 

For sufficiently large and fixed complex shifts, one distinctly perceives the number of iterations 
to gradually rise to a maximum around cr = — 4//i^. Upon reaching this maximum the number of 
iterations decreases slowly, exhibiting a steep descend for values of a around —8/h'^. The same 
observation has been reported in [22], where an eigenvalue analysis was performed to rigourously 
anticipate the Krylov convergence behaviour. Note that we have used standard Dirichlet boundary 
conditions, as opposed to the absorbing boundary conditions used in [22]. 

An explanation for the rise-and-fall of the iteration count can be found in the foundations of 
a general Krylov method, which typically traces and eliminates Ritz-eigenvalues starting from 
the most extreme eigenvalue in the spectrum. For wavenumbers around a = or a = -8//i^, the 
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problematic eigenvalues, i.e. the eigenvalues closest to the origin, are indeed extreme values as 
illustrated by Figure 14(a), implying the Krylov method will rapidly eliminate the problematic Ritz- 
eigenvalues which results in fast convergence. Figure 14(b) however shows that for cr around — 4//i^, 
the problematic eigenvalues are far from extreme, causing them to be eliminated very slowly by the 
Krylov method. This implies poor Krylov convergence, and many Krylov iterations are necessary to 
reach a sufficiently high precision on the solution. Additionally, since the complex Shifted Laplacian 
preconditioner shifts the spectrum over a fixed distance, preconditioning hardly improves the slow 
Krylov convergence. For additional information and a thorough analysis on the subject, we cordially 
refer the reader to [22]. 



4. CONCLUSIONS AND DISCUSSION 

In this paper we have analyzed the convergence of a complex Shifted Laplacian preconditioned 
Krylov solver for the Helmholtz problem. A Multigrid method is used to solve the Shifted Laplacian 
preconditioning system, which is a Helmholtz problem where the wavenumber /c^ is scaled by 
(1 + Pi). This results in an operator — A + a, where a = — + Pi) is complex valued. 

The asymptotic multigrid convergence rate of a two-, three- and four-grid scheme were analyzed 
theoretically with the aid of LFA. It is found that the convergence rate is mainly determined by 
the amplification factor maximum which appears at a single resonance frequency. The resulting 
convergence rate depends on the grid distance h, the complex shift parameter p and the wavenumber 
k. By increasing the complex shift /3, the maximum at the resonance decreases and the convergence 
rate improves. In general, the larger /3, the better the multigrid convergence rate. 

However, the larger we choose the value of /3, the further the complex shifted problem deviates 
from the original Helmholtz problem and the worse it will perform as a preconditioner. Indeed, a 
balance needs to be found between fast Krylov and multigrid convergence. 

From the expression of the convergence rate it is possible to define /3min as the smallest value of 
P for which the multilevel solution method of the shifted problem is stable. If P is taken smaller 
than /3min, unstable modes are bound to destroy the multigrid solver. 

When solving the preconditioner problem exactly using multigrid (up to discretization error 
order), the complex shift should always be taken larger than /3min- This ensures multigrid will 
converge when applied to the shifted problem. Experiments show that in this case the choice of 
P = Pmin is optimal, leading to a minimal number of Krylov iterations. However, when solving 
the preconditioner problem approximately using only a limited number of V-cycles, experimental 
results show that a minimum number of Krylov iterations is reached when choosing the shift 
parameter (3 slightly smaller than pmin- We conclude that /3min is a safe and near-optimal choice for 
the complex shift parameter, ensuring multigrid stability and solving the problem using a (nearly) 
minimal number of Krylov iterations. 

Additionally, we have shown that Prnin depends in an irregular way on the wavenumber k and 
the grid distance h, and that it is furthermore dependent on the number of pre- and post-smoothing 
steps. Choosing (3 around 0.5, as is common practice in the literature, ensures that the standard 
V(l,l)-cycle convergences for all kh < 0.625, since this (3 is distinctly larger than all corresponding 
Pmin- This indeed legitimizes the choice of (3 = 0.5 for practical purposes. 

For problems with regional space-dependent wavenumbers k{x), each region can be associated 
with a certain minimal complex shift Pmin- In that case, an indisputably safe choice for the complex 
shift parameter (3 would be the largest possible regional /3min appearing in the problem. 

Note that we have used homogeneous Dirichlet boundary conditions throughout the paper, while 
many applications use absorbing boundary conditions. This shifts the eigenvalues of the original 
Helmholtz problem away from the real axis and generally leads to better Krylov convergence. 
Consequently, the analysis in this paper is a discussion of the worst-case convergence scenario, and 
in practice much better Krylov convergence will be found than shown by the presented experiments. 
As Local Fourier Analysis makes no restictions on the boundaries however, the theoretical results 
presented within this text are generally valid. 
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